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O . Abstract. Using exact diagonalization calculations, we investigate the ground-state 

phase diagram of the hard-core Bose-Hubbard-Haldane model on the honeycomb 
lattice. This allows us to probe the stability of the Bose-metal phase proposed in 
Varney et ah, Phys. Rev. Lett. 107, 077201 (2011), against various changes in the 
originally studied Hamiltonian. 
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1. Introduction 



^ ■ In nature we are surrounded with examples of ordered phases at low temperatures — 

e.g. crystalline solid structures, magnetically ordered materials, superfluid and 
superconducting states, etc. While it is straightforward to think of these ordered phases 
melting as the temperature is increased into the familiar classical liquid or gaseous states 
that are commonplace in every aspect of our lives, it has been a long-standing question 
as to whether 'quantum melting' at zero temperature can act similarly to thermal effects 
and prevent ordering. For a quantum spin or boson system, the resulting state of matter 
is known as a quantum spin liquid [1]. The interest in such a hypothetical spin liquid has 
remained strong for decades, most prominently due to the discovery of high temperature 
superconductivity [2, 3]. 

Of critical importance is whether a two (or higher)-dimensional system can host a 
quantum spin liquid. At present, there exists a complete classification of quantum orders 
[4] , which divides hypothetical spin liquids into several distinct classes. Some theoretical 
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stability arguments have also been presented showing that there is no fundamental 
obstacle to the existence of quantum spin liquids [5]. Gapped spin liquid phases have 
been observed in dimer models [6-8], and also a family of special exactly-solvable toy 
models were discovered which can support gapped and gapless spin liquid phases [9]. 
Although these discoveries clearly demonstrated that a spin-liquid phase may appear 
in two (or higher) dimensions, at least in toy models, whether the same type of exotic 
phase can appear in a realistic spin system remains unclear. 

Very recently, there has been much numerical [10-17] and experimental [18-20] 
evidence to suggest the existence of gapped spin liquids in models with SU (2) symmetry, 
but it is still unclear why these simple models can support such exotic phases. Of 
particular note are the numerical discoveries of a gapped spin liquid in the Heisenberg 
model on the kagome lattice [11] and in the Hubbard model on a honeycomb lattice [10]. 
The existence of the latter is especially surprising and remains under debate [21]. The 
nature of this phase has been the subject of many works and it has been argued that 
next-nearest-neighbor exchange coupling is the mechanism responsible for the quantum 
spin liquid [13, 14, 22, 23]. However, despite a number of numerical investigations into 
this J1-J2 model, there are still open debates on whether the non-magnetic state present 
in this model is a valence bond solid [24-28] or a quantum spin liquid [12, 15, 16]. 

Gapless spin liquids, which may have low-lying fermionic spinon excitations that 
strongly resemble a Fermi-liquid state, have remained more elusive. Because of these 
excitations and because spin-1/2 models can be mapped onto hard-core boson models, 
some of these gapless spin liquids are often referred to as a Bose metal or Bose liquid. 
The hallmark feature of a Bose metal is the presence of a singularity in momentum 
space, known as a Bose surface [29-34]. However, unlike a Fermi liquid, where the 
Fermi wave vector depends solely on the density of the fermions, the Bose wave vector 
depends on the control parameters of the Hamiltonian and can vary continuously at 
fixed particle density. 

In this paper, we follow up on the proposition of such a putative Bose metal phase in 
a simple hard-core boson (XY) model on the honeycomb lattice [34], with an analysis of 
the stability of this phase against various changes in the Hamiltonian studied originally. 
First, we examine the dependence of the Bose wave vector on a (phase) parameter that 
makes the model transition between frustrated and non frustrated regimes. Next, we 
show that the phase identified as a Bose metal is stable to breaking of time-reversal 
symmetry and is present in the phase diagram of the hard-core Bose-Hubbard-Haldane 
(BHH) model, which features (at least) three phase transitions. 

The remainder of this paper is structured as follows. In section 2, we define the 
model Hamiltonian, briefly discuss the Lanczos algorithm, and define the key observables 
used in this study: the charge-density wave structure factor, the ground state fidelity 
metric, and the condensate fraction. Next, in section 3, we discuss the identifying 
characteristics of the Bose metal phase in the context of the hard-core boson (XY) model 
and show how the Bose wave vector evolves as the parameters are varied. In section 4, we 
discuss the three phase transitions that we can identify in the Bose-Hubbard-Haldane 
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model: BEC-CDW, BEC-Bose metal, and Bose metal(other phase)-CDW. The main 
results are summarized in section 5. 

2. Model and Methods 

2.1. Models 

The model proposed in [34] to exhibit a Bose metal phase is the spin- 1/2 frustrated 
antiferromagnetic-XF model on the honeycomb lattice 

H = Ji T,( S t$i + H.c.) + J 2 ^(S+Sj + H.c.) . (1) 

(ij) «W» 

where Sf is an operator that flips a spin on site % and J\ (J2) is the nearest- 
neighbor (next-nearest-neighbor) spin exchange. In this model, the next-nearest- 
neighbor coupling introduces frustration as long as J 2 > (antiferromagnetism). 

The Hamiltonian in (1) maps to a hard-core boson model (Sf — > b\, S~ — > b { , and 

Ji ->> U), 

H = t 1 J2(blbj + H.c.) + t 2 J2 (blbj + H.c). (2) 

Here b\ (b { ) is an operator that creates (annihilates) a hard-core boson on site i 
and ti (£2) is the nearest-neighbor (next-nearest-neighbor) hopping amplitude. This 
Hamiltonian was shown to feature four phases: a simple Bose-Einstein condensate 
(BEC) [a zero momentum (k = 0) condensate], a Bose metal (a gapless spin liquid), 
and two fragmented BEC states. The Bose metal (BM) was found to be the ground 
state of this model over the parameter range 0.210(8) < t 2 /ti < 0.356(9) [34]. 

To better understand the stability of the latter phase, we consider a strongly 
interacting variant of the Haldane model [35], the hard-core Bose-Hubbard-Haldane 
Hamiltonian [36] 

H = ti Y,( b l b j + H - c -) + *2 Y. {b > b :,<"'' + H - c + VY^n l n j , (3) 

(U> iv)) {ij) 

which reduces to (2) for 0jj = and V = 0. Here, V describes a nearest-neighbor 
repulsion and the next-nearest neighbor hopping term has a complex phase <pij = ±0, 
which is positive for particles moving in the counter-clockwise direction around a 
honeycomb. Note that the Hamiltonian in (3) can be mapped to a modified XXZ- 
model (Si~ -> b\, S~ -> 6 i: ni -> Sf + 1/2, t\ -> Ji, t 2 -> J2, and V -> J z ) 

H = J 1 Y.(S+Si + H.c.) + J 2 J2 {SfSje^ + H.c.) + J z £ (s* + ~) + ^ . (4) 

(ij) m (v) v 7 v 7 

The complex phase plays two important roles. Firstly, for <p ^ nir, time-reversal 
symmetry is explicitly broken. Therefore, we can use this control parameter to study the 
stability of the Bose metal phase against time-reversal symmetry breaking. Secondly, 
in the spin language, as we increase the value of <fi from to it, the sign for the 
next-nearest-neighbor spin-spin interaction is flipped from positive (<fr = 0) to negative 
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(0 = 7r), i.e., the next-nearest-neighbor spin exchange changes from antiferromagnetic 
to ferromagnetic. Since frustration in this model originates from the antiferromagnetic 
next-nearest-neighbor spin exchange, we can use to tune the system from a frustrated 
(0 = 0) to a non-frustrated (0 = it) regime, and thus it enables us to explore the role 
of frustration in stabilizing the BM phase. 

In what follows, t\ — 1 sets our unit of energy, and we fix t 2 = 0.3 to focus on 
transitions from the phase identified in [34] as a BM phase. This model has two limiting 
cases: (1) for V — > oo, the Ising regime, the ground state is a charge density wave 
(CDW) and (2) for V = and = ir, the non-frustrated regime, the ground state is a 
simple zero- momentum BEC with non-zero superfluid density (SF). 



2.2. Method and Measurements 

To determine the properties of the ground state of (3), we utilize a variant of the 
Lanczos method [37]. This technique provides a simple and unbiased way to determine 
the exact ground-state wave-function for interacting Hamiltonians. One limitation of 
the original algorithm is that the Lanczos vectors may lose orthogonality, resulting in 
spurious eigenvalues [38]. Orthogonality can be restored through reorthogonalization 
[39], which requires storing the Lanczos vectors. The storage needs can then be reduced 
utilizing a restarting algorithm, and the most successful techniques are the implicitly 
restarted [40, 41] and the thick- restart Lanczos algorithms [42]. These two methods 
are equivalent for Hermitian eigenvalue problems, and here we utilize the thick-restart 
method for its simplicity in implementation. 

A generic and unbiased way of determining the location of a quantum phase 
transition is related to the ground-state fidelity metric, g [36, 43-46]. The fidelity 
metric is a dimensionless, intensive quantity and is defined as 

_ 2 1-F(X,5X) 

where N is the number of sites and the fidelity F(X, 5X) is 

F(A,5A) = (v|>o(A)|*o(A + 5A)>, (6) 

where |\l/o(A)) is the ground state of H(X), and A is the control parameter of the 
Hamiltonian. 

For strong repulsive interactions, the ground state of the BHH model is a charge- 
density-wave (CDW) insulator, where one of the two sublattices is occupied while the 
other one is empty. This state spontaneously breaks the six-fold rotational symmetry 
down to three- fold but leaves the lattice translational symmetry intact. In addition, 
because of the diagonal character of the order established, the structure factor that 
describes this phase is maximal at zero momentum. Thus, we define the CDW structure 
factor Scow as 

^cdw = ^E(«-^)K-^)) ; (7) 

1,3 
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Figure 1: Clusters used in this study. 



where n® and nf are the number operators on sublattice a and b in the ith unit cell, 
respectively. 

Another possible ordered state is a Bose-Einstein condensate, where, in our 
model, bosons can condense into quantum states in which different momenta are 
populated. According to the Penrose-Onsager criterion [47], the condensate fraction 
can be computed by diagonalizing the one-particle density matrix pjj = (b\bj), 

fc = Ai/N b , (8) 

where Ai is the largest eigenvalue of pij and A& is the total number of bosons. In a BEC, 
the condensate occupation scales with the total number of bosons as the system size is 
increased, which is equivalent to stating that p^- exhibits off-diagonal long-range order 
[48]. Consequently, in a simple BEC, Ai ~ O(Nb) while all other eigenvalues are 0(1) 
[49]. Aside from a simple BEC, the eigenspectrum of the single-particle density matrix 
can signal fragmentation, where condensation occurs to more than one effective one- 
particle state [49, 50], and the Bose metal phase. In the former case, some of the largest 
eigenvalues are 0(iV&) and could even be degenerate. For the Bose metal, however, all 
of the eigenvalues of pij are ~ 0(1)- Thus finite size scaling of f c can help pinpoint the 
presence or absence of condensation. 

Further understanding of the latter two phases can be gained by calculating the 
single-particle occupation at different momentum points 

n(k) = (4« k ) + (4/? k ) , (9) 

where a k = J2ieA e * k Vi b\bi and /3 k = J2ieB e * k Vl b\bi are boson annihilation operators at 
momentum k for the A and B sublattices, respectively. In order to minimize finite-size 
effects and fully probe the Brillouin zone we average over 40 x 40 twisted boundary 
conditions [51, 52], 

<n(k)) fcA =fd0.f d6 y (n(k, 9 X , 9 y )) , (10) 
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Figure 2: (Color online) (a)-(b) Momentum distribution n(k) versus k [53] for the Bose 
metal phase in the XY model for tijt\ = 0.28 and 0.33, respectively. In both panels, 
40 x 40 twisted boundary conditions were averaged to generate n(k), and the Bose 
surface is indicated by a dashed red line, (c) The magnitude of the Bose-surface qs as 
a function of £2/^1 ■ 

where 9 a is the flux associated with the twisted boundary condition. 

For any finite-size calculation, there are a large number of clusters that one could 
study, each with slightly different symmetry properties. In this work, we focus solely 
on clusters that can be described by a parallelogram or "tilted rectangle" . The clusters 
used in this study are illustrated in figure 1 and are discussed in more detail in [36] and 
[34]. 

3. XY Model 

In a previous study [34], we reported that the phase diagram of the XY model (1) on a 
honeycomb lattice has three quantum phase transitions separating four distinct phases. 
The four phases are: (i) a BEC k = V (antiferromagnetism), (ii) a Bose metal (spin 
liquid), (iii) a BEC at k = M (a collinear spin wave), and (iv) a BEC at k = K (120° 
order) . 

The key signature of a Bose metal is the absence of any order and a singularity 
in the momentum distribution n(k). In figure 2(a) and 2(b), we show n(k) for two 
values of £2/^1 that are typical for the Bose metal phase. For this phase, n(k) features a 
^Ai-dependent Bose surface, which, as a guide to the eye, we indicate by a dashed red 
line. In general, the Bose wave vector qg at which the maxima of n(k) occurs increases 
with increasing £2/^1, as shown in figure 2(c). We emphasize that the maxima in n(k) 
do not reflect Bose-Einstein condensation as they do not scale with system size. 

4. BHH Model 



In this work, we present evidence that, in the (<f), V) plane [see (3)], the Bose-Hubbard- 
Haldane model for £2/^1 = 0.3 exhibits (at least) three phases at half-filling. For strong 
coupling V, the ground state is a charge-density- wave (CDW), while (at least) two 
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Figure 3: (Color online) (a) Fidelity metric g, (b) scaled structure factor L~ 7 ^S'cdw, 
and (c) scaled condensate fraction L 1+V f c as a function of interaction strength for various 
clusters with < 
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possible ground states exist at weak-coupling. In the frustrated regime (0 ~ 0) at 
V — 0, the system favors a Bose-metal, while the unfrustrated regime (<f) ~ ir) favors 
a BEC. Consequently, we find that there are three types of transitions: (i) BEC-CDW, 
(ii) BEC-BM, and (iii) BM(other phase)-CDW. 

Let us first consider the BEC-CDW transition driven by V at constant 0. In 
figure 3, we show the properties of the system for = n. In panel (a), we show the 
fidelity metric versus V, which has a smooth peak that grows with system size, indicative 
of a second-order phase transition (which would be unconventional in this case in which 
the system transitions between two ordered states) or a weakly first order transition. If 
the former is true, the structure factor would scale according to the rule: 

L-^S CDW = f[(V-V c )L^], (11) 

where N is the number of sites, L = N 1 ^ 2 is the linear dimension, and 7 = z/(2 — n). 
Because of our small lattice sizes, we cannot pinpoint the exact nature of this transition. 
For example, using a scaling analysis based on the 3D Ising [54] and XY universality 
classes [55] yield very similar results. In figure 3(b), we show the CDW structure factor 
scaled in accordance with the 3D XY universality class, resulting in V c = 3.71(7). 

We can check the robustness of this result by considering the condensate fraction, 
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which scales [56-58] according to 

LVc = »[(^-K)i 1/ 1, (12) 

where y = (d + z — 2 + rj). This is illustrated in figure 3(c), resulting in V c = 3.73(3). 
This result is quite close to the one obtained using the structure factor. We stress, once 
again, that although this appears to be a second-order transition between two ordered 
states, finite-size limitations do not allow us to rule out the possibility of a weak first- 
order transition or the existence of a small intermediate phase separating the BEC and 
CDW states. 

Next, we examine the properties of the model as one transitions from the Bose 
metal to the BEC state. In figure 4, we show the same quantities as in figure 3 (this 
time versus 0) for V = 0. The fidelity metric is plotted in figure 4(a), and peaks at 
approximately <j) ~ 0.88. In figure 4(b), we show the structure factor, which does not 
scale with finite size in either phase. Figure 4(c) depicts the scaled condensate fraction, 
yielding C = 0.84 ± 0.14, consistent with the peak in the fidelity metric. (Note that 
the 18A cluster experiences a level crossing for <fi < <p c .) As in the previous case, we 
cannot make definite statements about the nature of the transition between the Bose 
metal and the BEC state, but our results are consistent with a second order or a weakly 
first order transition. 

It is now interesting to study how the Bose surface changes as <fi increases and 
one transitions between the Bose-metal and the BEC phase. In figure 5, we show the 
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Figure 5: (Color online) Momentum distribution function n(k) versus k [53] for V = 
and (a) = 0, (b) = vr/9, (c) = vr/4, and (d) = tt/3. 



momentum distribution function for four values of and fixed V = 0. As seen in 
figure 5(b), the Bose surface reduces in size as departs from zero and continues to 
shrink until condensation occurs at k = T [see Figs. 5(c) and 5(d)] for > C . 

A third phase transition is expected as V is increased from zero and the BM phase 
is destroyed to give rise to the large V CDW phase. We illustrate this regime in figure 6 
by plotting the fidelity metric (left panels) and CDW structure factor (right panels) for, 
(a) = 0, (b) 7r/12, and (c) ir/6. In the left panels, for all three values of 0, one can 
see a sort of two-peak structure in the fidelity metric. 
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Figure 6: (Color online) Fidelity metric g (left) and structure factor S'cdw (right) as a 
function of interaction strength for various clusters with (a) = 0, (b) = 7r/12, and 
(c) = 7r/6. 
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Figure 7: (Color online) Momentum distribution function n(k) versus k [53] with 
constant = for (a) V = 0, (b) V = 0.25 and (c) V = 1. 

The large value of g at V = may be an indicator of a transition away from the 
Bose metal at V = + . It is somewhat similar to the behavior of g in both the one- 
dimensional Hubbard model, where the Mott-metal-insulator phase transition occurs for 
the onsite repulsion U = + [44], and the two-dimensional hole-doped t-J model [45], 
where <i-wave superconductivity was seen to develop for a superconducting inducing 
perturbation with vanishing strength. Another possibility is that the Bose-Metal is 
stable for positive and small values of V, but a transition to another phase occurs when 
V < 0. The peak produced by such a transition would also explain the structure we see 
in g. We have also investigated this model with negative values of V, and found that 
large peaks are present in the fidelity metric for V < 0. The position of those peaks 
had a strong dependence on the cluster geometry. Hence, exactly what happens to the 
BM phase in the region V ~ is something that requires further studies, maybe with 
other techniques that allow access to larger system sizes and a better finite size scaling 
analysis. 

For all clusters and values of depicted in the left panels in figure 6, one can also 
see a clear peak in the fidelity metric for finite values of V . This feature indicates the 
onset of CDW order. The structure factor S'cdw, depicted in the right panels in figure 6, 
make apparent that for values of V beyond that peak, the CDW structure factor scales 
with system size. 

In figure 7, we illustrate how the momentum distribution function changes in the 
presence of interactions at fixed = 0. n(k) is shown for V = in panel (a) and 
as the interactions are increased in panels (b) and (c). In figure 7(b), one can see 
that the Bose-surface broadens as V increases and moves closer to k = T. Increasing 
the nearest-neighbor repulsion further, so that the system enters in the CDW phase 
[figure 7(c)], results in a momentum distribution function that peaked at k = T, albeit 
without condensation. Instead, the structure factor S(k) is sharply peaked at k = T. 

A summary of our calculations for different values of V and is presented in 
figure 8 as the phase diagram of the hard-core BHH model at half- filling with t\ = 1.0 
and t-z = 0.3. For > 7r/4, the boundary of the CDW phase was identified by the 
crossing points in the scaling of the structure factor (figure 3). The boundary of the 
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Figure 8: (Color online) Phase diagram for the Bose-Hubbard-Haldane model with 
parameters t\ = 1.0 and t 2 = 0.3. The solid red squares are determined by the crossing 
point in the scaling of Scdw- The solid blue circles are from the crossing point in f c . 
The green triangles are the average of the location of peak in the fidelity metric for 
the largest system sizes. The Bose Metal (BM) phase is indicated by the thick, dashed 
magenta line. 



BEC phase was identified by the crossing points in the scaling of the condensate fraction 
(Figs. 3 and 4), and, for small values of V, also using the maximum of the fidelity metric 
for the largest systems sizes (figure 4). For <p < 7r/4, the CDW transition boundary was 
determined by the position of the maximum in the peak in the fidelity metric for the 
largest system sizes (figure 6). Note that, in that regime, the Bose metal phase was 
found to be stable for V = 0. On the other hand, for V between and the boundary of 
the CDW phase, the large value of the fidelity metric, as well as the behavior of several 
observables studied in that region, prevent us from making a clear statement about the 
nature of the ground state. 

5. Conclusion 

In summary, we have studied the phase diagram of the hard-core Bose-Hubbard-Haldane 
Hamiltonian, which has allowed us to probe the effect of perturbations on the Bose metal 
phase found in the frustrated XY model on a honeycomb lattice [34]. In particular, we 
explored the parameter dependence of the Bose wave vector and verified that the Bose 
metal is stable under the effects of time-reversal and chiral symmetry breaking. We 
identified three phases in the phase diagram of the BHH model; (I) a Bose metal, 
(II) a BEC, and (III) a CDW. The phase transitions between the different phases 
were identified utilizing the ground state fidelity metric, the CDW structure factor, 
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the condensate fraction, and the momentum distribution. 

The BEC-CDW transition appears to be second order, although finite-size effects 
prevent us from ruling out the possibility of a weak first-order transition or the existence 
of an intermediate phase separating the BEC and the CDW states. If this transition 
is indeed a direct second order phase transition, the critical point would be highly 
nontrivial, and could be an example of deconfined criticality. 

We have also found that the Bose metal is destroyed upon increasing V, before 
the Heisenberg point for nearest neighbor interactions V = 2J\ can be reached. The 
presence of a next-nearest-neighbor repulsion may change this and result in transitions 
to other exotic phases. 
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